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Abstract. Knowledge of the structure of the coronal mag- 
netic field is important for our understanding of many solar 
activity phenomena, e.g. flares and CMEs. However, the 
direct measurement of coronal magnetic fields is not possi- 
ble with present methods, and therefore the coronal field has 
to be extrapolated from photospheric measurements. Due to 
the low plasma beta the coronal magnetic field can usually 
be assumed to be approximately force free, with electric cur- 
rents flowing along the magnetic field lines. There are both 
observational and theoretical reasons which suggest that at 
least prior to an eruption the coronal magnetic field is in a 
nonlinear force free state. Unfortunately the computation of 
nonlinear force free fields is way more difficult than poten- 
tial or Unear force free fields and analytic solutions are not 
generally available. We discuss several methods which have 
been proposed to compute nonlinear force free fields and fo- 
cus particularly on an optimization method which has been 
suggested recently. We compare the numerical performance 
of a newly developed numerical code based on the optimiza- 
tion method with the performance of another code based on 
an MHD relaxation method if both codes are applied to the 
reconstruction of a semi-analytic nonlinear force-free solu- 
tion. The optimization method has also been tested for cases 
where we add random noise to the perfect boundary con- 
ditions of the analytic solution, in this way mimicking the 
more reahstic case where the boundary conditions are given 
by vector magnetogram data. We find that the convergence 
properties of the optimization method are affected by adding 
noise to the boundary data and we discuss possibiUties to 
overcome this difficulty. 



1 Introduction 

The magnetic field is the most important quantity for under- 
standing the plasma structure of the solar corona and the ac- 
tivity phenomena it shows. Unfortunately a systematic direct 
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measurement of the coronal magnetic field is not feasible 
with the observational methods presently available. There- 
fore the extrapolation of the magnetic field into the corona 
from measurements of the magnetic field at the photospheric 
level is extremely important. 

It is generally assumed that the magnetic pressure in the 
corona is much higher than the plasma pressure (small plasma 
/3) and that therefore the magnetic field is nearly force-free 
(for a critical view of this assumption see Gary , 2001). The 
extrapolation methods based on this assumption include po- 
tential field extrapolation (e.g., Schmidt , 1964; Semel , 1967), 
linear force-free field extrapolation (e.g., Chiu and Hilton 
, 1977; Seehafer , 1978, 1982; Semel , 1988) and nonlin- 
ear force-free field extrapolation (e.g., Amari et al. , 1997). 
Methods for non-force-free field extrapolation have also been 
developed (Petrie and Neukirch , 2000) but are not used rou- 
tinely. 

Whereas potential and linear force-free fields can be used 
as a first step to model the general structure of magnetic 
fields in the solar corona, the use of non-linear force free 
fields is essential to understand eruptive phenomena as there 
are both observational and theoretical reasons which suggest 
that the pre-emptive magnetic fields are non-linear force-free 
fields (see below). The calculation of nonlinear force free 
fields is complicated by the intrinsic nonlinearity of the un- 
derlying mathematical problem. Another problem which be- 
comes especially important when nonlinear force free field 
are used for magnetic field extrapolation is the correct formu- 
lation of the problem with respect to boundary values. The 
available magnetograph observation provide either the line- 
of-sight magnetic field (Bios), which is sufficient for poten- 
tial and hnear force free fields, or all three components of 
the photospheric magnetic field. Only the latter information 
is sufficient to determine non-linear force free fields com- 
pletely. There are, however, tremendous difficulties to be 
overcome both in extracting the necessary information about 
the magnetic field components on the boundary from the data 
in an unambiguous way and in the use of this information in 
the computation of the corresponding magnetic field. Even 
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though a lot of research has been devoted to these problems 
in the past so far none of the proposed methods has been 
found to be outstanding in combining simpUcity, usability 
and robustness. 

The purpose of the present contribution is to assess some 
of the methods presently available for calculating non-linear 
force-free fields, with a view to implementing them into a 
general extrapolation code which also takes stereoscopic in- 
formation into account in the extrapolation process. A first 
version of such a method based on linear force free extrapola- 
tion has recently been proposed (Wiegelmann and Neukirch 
, 2002). As stated in Wiegelmann and Neukirch (2002) and 
bearing in mind the remarks above, it is would be very de- 
sirable to generalize the extrapolation method to nonlinear 
force free fields. 

As a first step in this direction we investigate the robust- 
ness of the optimization method recently proposed by Wheat- 
land et al. (2000), and compare its performance with another 
method based on MHD relaxation. The paper is organized as 
follows. In Sect. 2 we discuss the basic equations and back- 
ground theory. Section 3 describes various methods which 
have been proposed for nonlinear force free field extrapola- 
tion. Out of these we focus on the optimisation method and, 
for the case of analytically given boundary data, on the MHD 
relaxation method. These methods are tested using an analyt- 
ical solution in Sect. 4. The optimisation method is also sub- 
jected to tests with more general boundary conditions where 
we add noise to the given analytical boundary data. The pa- 
per closes with the conclusions in Sect. 5. 

2 Basic equations: Force-Free Equilibria 

For force free fields the equations to solve are 

jxB = (1) 
VxB = Moj, (2) 
V-B = 0. (3) 

Equation (1) implies that for force free fields the current den- 
sity and the magnetic field are parallel, i.e. 

Moj = aB. (4) 

Here a is a function of space. This function has to satisfy the 
equation 

B • Va = (5) 

which is obtained by taking the divergence of Eq. (4) and 
making use of Eqs. (2) and (3). Substituting Eq. (4) into Eq. 
(2) we get a second equation which determines B : 

V X B = aB (6) 

Equation (5) implies that a is constant along magnetic 
field lines, but that it can vary across the magnetic field. Ob- 
viously, Eq. (5) is solved by a = implying that j = by 
Eq. (4). In this case we obtain potential fields. Another ob- 
vious possibility is a constant, but nonzero. This is the case 



of linear (or constant-a) force free fields. Standard methods 
for calculating potential and linear force free field are avail- 
able (see e.g. Gary , 1989) and they can at least give a rough 
impression of the coronal magnetic field structure. 

There are, however, both observational and theoretical ar- 
guments that at least the magnetic field prior to eruptive pro- 
cesses in the corona is not a Unear force free (or potential) 
field. If the normal component or the line-of-sight compo- 
nent of the magnetic field on the boundary is given the cor- 
responding potential field is uniquely determined and is ac- 
tually the magnetic field with the lowest magnetic energy for 
this boundary condition. Since magnetic energy is needed 
to drive an eruptive process pre-emptive coronal magnetic 
fields cannot be expected to be in a potential state. This is 
also corroborated by observations (e.g. Hagyard , 1990; Hag- 
yard et al. , 1990; Falconer et al. , 2002). Similar observa- 
tions also indicate that it will usually be difficult to match the 
boundary magnetic field with a linear force-free field. Also, 
from a more general point of view it seems highly unhkely 
that the complex pre-eruption coronal magnetic field which 
is slowly but constantly stressed by changing boundary con- 
ditions will have the current density distribution of a hnear 
force-free field, i.e. a constant a. A (strongly) localized cur- 
rent density in the erupting configuration seems to be a lot 
more natural. Another argument which is sometimes put for- 
ward to rule out linear force-free magnetic fields for models 
of pre-emptive states is that linear force free fields minimize 
the magnetic energy under the assumption of global helicity 
conservation (e.g., Taylor , 1974), and therefore if helicity 
conservation could be assumed to hold for a coronal erup- 
tion, a hnear force-free field would not be able to provide the 
energy for the emption. Even though it has been proposed to 
apply variants of Taylor's hypothesis to the state of the coro- 
nal magnetic field on larger scales (e.g., Heyvaerts and Priest 
, 1984), it generally does not seem to apply for models of 
eruptive processes (e.g., Amari and Luciani , 2000). Further- 
more, there is the possibility of helicity transport into and out 
of the region of interest and for these reasons, we do not want 
to put too much emphasis on this argument. 

However, we beheve that the other arguments are strong 
enough to assume that generally one has to take into account 
that for the force free magnetic fields of pre-emptive states 
the value of a changes from field fine to field fine. This au- 
tomatically leads us to consider nonlinear force free fields. 

3 Extrapolation Methods 

The problem we discuss in this section is to compute non- 
linear force free magnetic fields in Cartesian coordinates if 
given all three components of the magnetic field on the bound- 
ary, i.e. the photosphere. Various methods have been pro- 
posed (see e.g. Amari et al. , 1997; Demoulin et al. , 1997; 
McClymont et al. , 1997; Semel , 1998; Amari et al. , 2000; 
Yan and Sakurai , 2000) of which we will discuss three. For 
the first method we assume that the magnetic field is given 
on the lower boundary (photosphere), whereas for the sec- 
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ond and third method we assume that B is prescribed on the 
surfaces of a volume V. 

3 . 1 Direct Extrapolation 

Direct extrapolation (Wu et al. , 1990) is a conceptually sim- 
ple method, in which the Eqs. (6) and (5) are reformulated 
in such way that they can be used to extrapolate the photo- 
spheric boundary conditions given by a vector magnetogram 
into the solar corona. 

If we suppose that the boundary condition on the lower 
boundary is given by Bo(a;, y, 0) then we can calculate the 
^-component of the current density at z = by using Eq. (2) 



dB 



yO 



dx 



dy 



(7) 



Knowing the z-component of the photospheric magnetic field 
(Bzo), we can use Eqs. (7) and (4) to determine a{x, y, 0) = 



J 

B 



(8) 



zO 



We remark that special care has to be taken at photospheric 
polarity inversion Unes, i.e. lines along which B^o = (see 
e.g., Cuperman et al. , 1991). Equation (8) allows us to cal- 
culate the X and y-component of the current density using 
Eq. (4) 



jyo 



ao Byo- 



(9) 
(10) 



We now use Eq. (3) and the x and y-component of Eq. (2) 
to obtain expressions for the ^-derivatives of all three mag- 
netic field components in the form 



dBxo 

dz 
dByo 

dz 
dB,o 
dz 



jyO + 
dB,o 

dy 

dBxo 
dx 



dB,o 
dx 



Jxa, 
dByo 



dy 



(11) 
(12) 

(13) 



The idea is to integrate this set of equations (numerically) in 
the direction of increasing z, basically repeating the previous 
steps at each height. 

Unfortunately it can be shown that the formulation of the 
force free equations in this way is unstable because it is an 
ill-posed problem (e.g., Cuperman et al. , 1990; Amari et al. 
, 1997). In particular one finds that exponential growth of the 
magnetic field with increasing height is a typical behaviour. 
The reason for this is that the method transports information 
only from the photosphere upwards. Other boundary condi- 
tions, e.g. at an upper boundary, either at a finite height or 
at infinity cannot be taken into account. Attempts have been 
made to regularize the method (e.g., Cuperman et al. , 1991; 
DemouUn et al. , 1992), but cannot be considered as fully 
successful. 



3.2 MHD-Relaxation 

Another possibility to calculate nonlinear force free fields is 
by the method of MHD relaxation (e.g., Chodura and Schliiter 
, 1981). The idea is to start with a suitable magnetic field 
which is not in equilibrium and to relax it into a force free 
state. This is done by using the MHD equations in the fol- 
lowing form : 

z/v = (V X B) X B (14) 

(15) 



E-F V X B = 
dB 

'dt 

V B = 



-V X E 




(16) 
(17) 



The equation of motion (14) has been modified in such a way 
that it ensures that the (artificial) velocity field is reduced. 
Equation (15) ensures that the magnetic connectivity remains 
unchanged during the relaxation. The relaxation coefficient 
can be chosen in such way that it accelerates the approach to 
the equiUbrium state. We use 

iy = -\Bf (18) 

with fj, = constant. Choosing the relaxation coefficient v pro- 
portional to B^ speeds up the relaxation process in regions 
of weak magnetic field (see Roumeliotis , 1996). Combin- 
ing Eqs. (14), (15), (16) and (18) we get an equation for the 
evolution of the magnetic field during the relaxation process : 

dB _ / [(VxB)xB]xB - 
^-mVx I 



(19) 



This equation is then solved numerically starting with a given 
initial condition for B. Equation (19) ensures that Eq. (17) 
is satisfied during the relaxation if the initial magnetic field 
satisfies it. 

The difficulty with the relaxation method in the present 
form is that it cannot be guaranteed that for the boundary 
conditions we impose and for a given initial magnetic field 
(i.e. given connectivity), a smooth force-free equilibrium ex- 
ists to which the system can relax. If such a smooth equi- 
librium does not exist the formation of current sheets is to be 
expected which will lead to numerical difficulties. Therefore, 
care has to be taken when choosing an initial magnetic field. 
This clearly limits the applicability of the relaxation method, 
and further work will be needed to overcome this obstacle. 
We therefore only apply the relaxation method for compar- 
ison with the optimization method to cases where we know 
the boundary conditions to be compatible. 

One can show, however, that the relaxation method con- 
verges under less restrictive boundary conditions. By multi- 
plying Eq. (16) by B, integrating over the complete compu- 
tational volume V and using Eq. (15), we find that 

^WB=(f [iv-B)B-B^v]-ndS- [ iixB)-vdV,(20) 

where 



dV 



Wb 



2X 



B^dV. 



(21) 
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Under line-tying boundary conditions we have v = on the 
boundary, and only the volume integral survives. Using Eq. 
(14) and the fact that u > one can show that in this case 



^ (j X BfdV < 0, 



(22) 



and an equiUbrium is achieved if and only if the magnetic 
field is force free. 

We emphasize that these boundary conditions are usually 
not compatible with prescribing all three magnetic field com- 
ponents on the boundaries, but only the normal component. 
Therefore, in the following we discuss the relaxation method 
only for test cases in which we know that the boundary con- 
ditions are compatible. This is done for comparison with the 
other method we discuss now, the optimization method. 

3.3 Optimization 

Recently Wheatland et al. (2000) have proposed an opti- 
mization method which refines a proposal by Roumehotis 
(1996). In this method the functional 



L = 



[ [B-2|(VxB) xBp + |V-Bp] d^y 
Jv 



(23) 



is minimized. Obviously, L is bounded from below by 0. 
This bound is attained if the magnetic field satisfies the force 
free equations 

(V X B) X B = (24) 

V-B = 0. (25) 

The variation of L with respect to an iteration parameter t 
leads to (see Wheatland et al. (2000) for details) 



dB 



dt 



■Fd^V 



dB 



dt 



■Gd^S 



(26) 



1 dL _ 
2~dt ~ 
where 

F = V X (f2 X B) - f2 X (V X B) 

-v(n ■ B) -I- n(v B) + n^B (27) 

G = n X (12 X B) - n{n ■ B) (28) 

ft = B-^ [(V X B) X B - (V • B) B] (29) 

The surface term in (26) vanishes for dB/dt = ' on the 
boundary and L decreases for 



(30) 



This leads to an iteration procedure for the magnetic field 
which is based on the equation 



dB 

'dt 



[(V X B) X B] X B 



-I- n |-f2 X (V X B) - V(f2 • B) 
-|-n(V ■ B) -I- fi^ b| 



(31) 



^This condition makes it necessary that all three components of the mag- 
netic field have to be prescribed on the six bovmdaries of the computational 
box. In section 4.2.3 we discuss how the surface integral in (26) can be used 
to update the lateral and top boundary conditions during the optimization 
process. Under this condition only the bottom boundary is prescribed dme 
independently with the help of photospheric vector magnetograms. 



Here is a constant which can be chosen to speed up the 
convergence of the iteration process. Equation (19) and the 
leading terms of Eq. (31) are identical, but Eq. (31) contains 
additional terms. 

For this method the vector field B is not necessarily sole- 
noidal during the computation, but will be divergence-free if 
the optimal state with L = is reached. A disadvantage of 
the method is that it cannot be guaranteed that this optimal 
state is indeed reached for a given initial field and boundary 
conditions. If this is not the case then the resulting B will 
either be not force free or not solenoidal or both. 



4 Tests 

We encoded the MHD-relaxation method and the optimiza- 
tion method in one program. The code is written in C and 
uses 4th order finite differences on an equal spaced grid. The 
time iteration is computed with the method of steepest gradi- 
ent (see e.g. text books like Geiger and Kanzow , 1999). The 
program is parallehzed with OpenMR 

4. 1 The semi-analytical test field 

To test the reconstruction methods, we try to reconstruct a 
semi-analytic nonlinear force free solution found by Low and 
Lou (1990). Wheatland etal. (2000) have used similar tests. 
The main difference between their paper and ours is in the 
diagnostic quantities used and in the comparison with i) the 
relaxation method and ii) the noisy boundary data (which we 
think is more representative of a realistic situation). Low and 
Lou (1990) solved the Grad-Shafranov equation for axis- 
symmetric force free fields in spherical coordinates r, 6, (p. 
For axis- symmetry the magnetic field can be written in the 
form 

1 _ . ^^2) 



B = 



r smf 



1 dA dA 
r do dr 



where A is the flux function, and Q represents the t/i-component 
of B, depending only on A. The flux function A satisfies the 
Grad-Shafranov equation 



d^A ^ l-ii^ d^A ^ ^dQ _^ 

gj.2 J.2 Q^2 



(33) 



where /x = cos 6. Among others. Low and Lou (1990) derive 
solutions for 



= a = constant 



dQ 
dA 

by looking for separable solutions of the form 



(34) 



(35) 



The solutions are axis- symmetric in spherical coordinates 
with a point source at the origin, but if used for testing force 
free codes in Cartesian geometry the symmetry is no longer 
obvious after a translation which places the point source out- 
side the computational domain and a rotation of the syname- 
try axis with respect to the Cartesian coordinate axis. 
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4.2 Test rans 

For the test we have first calculated Bz{x, y, 0) on the pho- 
tospheric boundary z = for the appropriate Low and Lou 
solution. Using this as a boundary condition we calculate 
the potential field corresponding to this boundary conditions 
inside our computational box. The potential field is com- 
puted with help of a method developed by Seehafer (1978). 
This method gives the components of the magnetic field in 
terms of a Fourier series. The observed magnetogram (or 
here the extracted magnetogram from the Low and Lou so- 
lution) covers a rectangular region extending from to Lj; 
in X and to Ly in y is artificially extended onto a rectan- 
gular region covering —L^ to and —Ly to Ly by taking 
an antisymmetric mirror image of the original magnetogram 
in the extended region, i.e. Bz{—x,y) = —Bz{x,y) and 
Bz{x,—y) = —Bz{x,y). We use a Fast Fourier Transfor- 
mation (FFT) scheme (see Alissandrakis (1981)) to deter- 
mine the coefficients of the Fourier series. For more details 
regarding this method see Seehafer (1978, 1982). 

The potential field is used as starting field inside the com- 
putational box but the Low and Lou B field is imposed on 
the boundaries. We then use either the MHD relaxation or 
the optimization method to calculate the correct force free 
equilibrium field. During the computations we calculate the 
quantities L/[T^m] (for both relaxation and the optimiza- 
tion method), the absolute value of the Lorentz force |J x 
'Q\/[nN m~^] (averaged over the numerical grid), the value 
of |V-B|/[mGaussm~^] (averaged over the numerical grid), 
and the difference between the numerical magnetic field and 



the known analytical solution 



(averaged over 



|B(t)-B„„„|-' 

the numerical grid) at each time step. In section 4.2.1 and 
4.2.2 all components of the magnetic field are fixed on all six 
boundaries. In section 4.2.3 only the bottom boundary con- 
dition is prescribed time-independent and the lateral and top 
boundaries are updated during the iteration. The details are 
described in 4.2.3. 

4.2.1 Standard tests 

In Fig. 1 we show three-dimensional plots of selected mag- 
netic field lines for the Low and Lou solution, the potential 
field calculated by taking the Low and Lou B^ on the lower 
boundary, the field of the MHD relaxation method after 50, 
500 and 5. 000 relaxation time steps, and the same plots for 
the optimization method. In these runs we used a grid size 
of 40 X 40 X 20. The colour coding of the bottom boundary 
indicates the Bz distribution on that boundary. 

It can bee seen that the potential field which is used as 
starting field of the computations for both methods is clearly 
different from the Low and Lou test field. The state of the 
system after 50 steps still shows some resemblance to the 
initial potential field for both methods but the field Unes have 
started to evolve away from the potential field. One can 
also see small differences between the two methods. Af- 
ter 5, 000 steps no obvious differences between the magnetic 
field reached with either method and the Low and Lou field 



Table 1. Details of runs to reconstruct a Low and Lou (1990) solution. 
The values of the parameters used by Low and Lou are / = 0.3 and $ 

■ The first column contains the used grid size and comments which code 
(relaxation or optimization) and boundary conditions have been used. If 
not specified the boundary conditions have been extracted from the analytic 
solution. We specify if noise is added to the boundary conditions. The 
remark Pot-boundary means that the lateral and top boundary conditions are 
not extracted from the analytic solution, but prescribed as a global potential 
field. The second column contains the iteration step. The third column 
contains the value of the functional L, the fourth column the Lorentz force 
(averaged over the numerical grid) and the last column the relative error 
compared with the analytic solution. 



rix X Uy X Uz 


Step 


L 




Relative Error 






40 X 40 X 20 


discret. error 


3.2 


2.0 


Reference 


Start 





46573 


264 


0.42 


Relaxation 


50 


3092 


38.2 


0.26 




500 


769 


7.5 


0.095 




5,000 


25.0 


1.4 


0.0016 


Optimization 


50 


832 


44.1 


0.19 




500 


23.1 


7.0 


0.03 




5,000 


3.15 


2.05 


0.0001 


1% noise 


5,000 


9.5 


3.7 


0.00017 


10% noise 


5,000 


645.0 


29.2 


0.0055 


20% noise 


5,000 


2584 


58.6 


0.028 


Pot-boundary 


500 


468 


9.1 


0.39 


Pot-boundary 


5,000 


458 


6.3 


0.37 


80 X 80 X 40 


discret. error 


0.6 


0.65 


Reference 


Start 





109203 


297 


0.48 


Optimization 


50 


1590 


52 


0.26 




500 


81 


12 


0.079 




5,000 


0.79 


1.02 


0.0013 




25,000 


0.6 


0.65 


0.000003 


10% noise 


25,000 


1876 


37.6 


0.0078 


Pot-boundary 


1,000 


563 


12.8 


0.42 


Pot-boundary 


10,000 


529 


5.3 


0.36 



can be seen. 

To quantify this statement we show in Fig. 2 a comparison 
of the evolution of the four diagnostic quantities for the two 
methods during the computation. For both methods all four 
diagnostic quantities decrease during the computation but af- 
ter 5, 000 steps the optimization method shows significant 
smaller values for the quantities L, V • B and the compari- 
son with the analytic solution. Therefore we can state that in 
this case the optimization method seems to be slightly more 
efficient than the MHD relaxation method. This is corrobo- 
rated by a look a Table 1 in which we summarize the main 
results of the various test runs we have made. The first row 
shows the discretisation errors for L and the Lorentz force 
if the known solution is discretised on a 40 x 40 x 20-grid 
and used to calculate the values of this quantities. The sec- 
ond row contains the values of L, the Lorentz force and the 
relative error after the interior grid points have been replaced 
by the potential field calculated from the photospheric B^ of 
the Low and Lou solution. The relatively large values of the 
three diagnostic quantities show the deviation from the equi- 
librium. 

The next two rows show how the three diagnostic quanti- 
ties evolve during the MHD relaxation method for the 40 x 
40 X 20-grid. One can see that after 5, 000 steps the value 
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Fig. 1. Top row, left panel: A set of selected field lines of the Low and Lou solution. Top row, right panel: The same field lines for the corresponding potential 
field calculated as described in the main text. The difference between the two fields is obvious. Middle row: The field lines after 50 (left) and 5,000 (right) 
steps of the MHD relaxation method. Bottom row: The field lines after 50 (left) and 5,000 (right) steps of the optimization method. The box drawn shows 
the spatial extension of the numerical box. The colour coding at the bottom of the boxes shows the normal component of the photospheric magnetic field in 
Gauss. 
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51 ' ' 

4 h 

3 h-. ...... 




2000 3000 4000 



Steps 



Log (iJxBl) 

2.5 F ' ' ' 




1000 2000 3000 4000 

Steps 



Loc(l7-H 




-fiO F 

1000 2000 3000 4000 

Steps 




1000 2000 3000 4000 

Steps 



Fig. 2. Evolution of the diagnostic quantities. Top panel: functional L 
(defined in Eq. (23)), second panel: Lorentz force | j X B | (averaged over 
the numerical grid), third panel: | V • B| (averaged over the numerical grid) 

bottom panel: relative error J — — "p"' (averaged over the numerical 
grid). The quantities are drawn for the MHD relaxation method (dotted 
lines) and for the optimization method (solid lines). The grid size used was 
40 X 40 X 20. 



of L has dropped by almost three orders of magnitude, but is 
still one order of magnitude above the value calculated with 
the discretised exact solution. The Lorentz force has dropped 
to the level of the dicretisation error, and the relative error has 
dropped more than two orders of magnitude. 

We are giving the same quantities for the optimization 
method in the following three rows. It can be seen that espe- 
cially the values of L are always way below the correspond- 
ing values of the MHD relaxation method. This is not sur- 
prising as the optimization method relies on the minimization 
of the functional L. The final value of L is actually below the 
discretisation error. It can also be seen that the relative error 
after 5, 000 iterations is more than one order of magnitude 
smaller than the corresponding value of the MHD relaxation 
method. 

We have made a step towards checking numerical conver- 
gence for the optimization method by repeating the calcula- 
tion on a grid with doubled resolution (80 x 80 x 40). The 
results are shown in the lower part of Table 1 . The first thing 
to notice is that the discretisation error is almost an order of 
magnitude smaller than for the previous grid. For this grid 
25, 000 iteration steps have been carried out, and after that 
the values of L and the Lorentz force have reached the level 
of the discretisation error. The relative error is more than two 
orders of magnitude smaller than for the coarser grid. 

4.2.2 Effect of adding noise 

The previous calculations have been carried out under the 
assumption that the magnetic field on the boundary of the 
computational box is known exactly. Such an ideahzed situ- 
ation will not be found when real data are used. Vector mag- 
netograms will have finite resolution and suffer from obser- 
vational uncertainties making the reconstruction potentially 
more difficult. Especially the optimization method has been 
proposed in view of coping with these difficulties (Wheat- 
land et al. , 2000). 

To investigate how the optimization method works if the 
boundary conditions are not given by the exact analytical 
field, but to keep control over the amount of uncertainty, 
we have carried out test runs with the optimization method 
adding random noise with to the boundary conditions. We 
add the noise by multiplying the exact boundary conditions 
with a number ^ where ^ is a random number in the range 
—ni < S < ni and ni is the noise level. 

To study the effect of the noise we have done runs with 
different ampUtudes of noise on the 40 x 40 x 20 grid. The 
evolution of L and the Lorentz force with the numbers of it- 
eration for various noise levels are shown in Fig. 3. It can 
clearly be seen that the method converges less and less well 
with increasing noise. For noise levels of 10 % and 20 %, the 
corresponding values after 5, 000 steps for the 40 x 40 x 20 
grid are given in Table 1. We have also carried out a run on 
the 80 X 80 X 40 grid with a noise level of 10 %. For this run 
the values of the diagnostic quantities after 25, 000 steps are 
still higher than the values of the corresponding quantities on 
the coarser grid after 5, 000 steps. We have to conclude that 
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Fig. 3. Evolution of L and |j x B| for the optimization code for various 
levels of boundary value noise. The lowest Unes correspond to an ideal 
vector magnetogram and the higher lines to vector magnetograms with 1%, 
10% and 20% noise, respectively. The grid size used here is 40 x 40 x 20. 
The convergence of the method gets smaller with increasing noise level. 

even a 10 % uncertainty in the boundary conditions could af- 
fect the convergence of the optimization method quite badly. 
The main problem is that with noise and fixed boundary con- 
ditions on all six boundaries the boundary condition are over- 
imposed and will no longer be compatible. We discuss how 
this problem can be solved in the next section. 

4.2.3 The problem of the lateral and top boundary condi- 
tions 

Until now the runs have been carried out under the assump- 
tion that the magnetic field boundary conditions are well known 
on all six boundaries of the box. Unfortunately real vector- 
magnetograms only provide information regarding the pho- 
tospheric magnetic field. The top and side boundary con- 
ditions are unknown and have to be fixed somehow. Here 
we want to investigate the influence of the choice of these 
boundary conditions. An additional problem occurs by fix- 
ing the vector magnetic field on all six boundaries (as done in 
the previous sections) because one overimposes the bound- 
ary conditions. For the test cases in the previous sections 
the boundary conditions have been extracted from an ana- 
lytic solution and consequently the boundaries are automati- 
cally compatible (this implies that there are no problems with 
over-imposed boundary conditions). Section 4.2.2 showed 
that this compatibility might get lost if the boundaries con- 
tain noise. To overcome this difficulty and to take into ac- 
count that measured data are only available for the bottom 
boundary we describe how the strict condition ^ = on 
all boundaries can be avoided within the optimization proce- 
dure. First we have to impose the boundary conditions as a 



well posed problem. A popular weU posed boundary condi- 
tion for non linear force free fields is to impose the normal 
component i?„ of the magnetic field and the normal com- 
ponent of the current density j,, in regions for a positive (or 
negative) Bn- (see Aly , 1989, regarding the compatibility of 
photospheric vector magnetograph data). These conditions 
have been derived under the very strict constraint of a flux 
balanced magnetogram where all magnetic flux is closed. 
This implies that each magnetic field line starting at one point 
on the photosphere also has to end on the photosphere (in a 
region of opposite magnetic flux, but with the same value of 
a). Such strictly isolated active regions seem to be a too re- 
strictive constraint for real vector magnetograms. Real vec- 
tormagnetograms provide all three components of the mag- 
netic field for either sign of -B„ and thus also j„ on the com- 
plete bottom boundary. It seems reasonable to impose these 
observed data on the bottom boundary. The price we have 
to pay is that field lines starting on the photosphere might 
pass the lateral and top boundary. ^ Consequently we need 
to update the lateral and top boundary during the iteration. 

As a first step towards consistent boundary conditions for 
the optimization method, we choose a potential field on the 
lateral and top boundaries. Firstly, a potential field can be 
easily computed just from the measured normal component 
of the photospheric magnetic field and secondly, it can be 
assumed that the solar magnetic field is reasonably well ap- 
proximated by a potential field outside active regions. Here 
we compute the potential field on the boundary as a global 
potential field computed from the distribution at ^; = 
alone. Table 1 shows that the value of L for the 40 x 40 x 20 
grid drops to a slightly lower value than for the correspond- 
ing run with 10% noise, while the remaining |j x B| forces 
are a factor of three lower. The relative error compared with 
the exact solution is quite high here, which is no wonder be- 
cause the magnetic field is forced to stay potential close to the 
side and top boundaries. One might notice that the error in 
the forces is a factor of 3 larger than the discretisation error, 
while L is more than two orders of magnitude larger than the 
discretisation error. The reason is that an inconsistency oc- 
curs at the edges of the box between the photospheric bound- 
ary and the side boundaries. The Low and Lou solution is not 
potential on the photosphere close to the side boundaries, but 
the chosen side boundary conditions are potential. 

To overcome this difficulty we cannot fix the lateral and 
top boundaries during the iteration, but have to relax also 
these boundaries. As the iteration equation (30) has been 
derived under the condition ^ = on all boundaries we 
have to extend the optimization approach. If we allow ^ ^ 
on some boundaries the surface term in equation (26) does 



^If the flux is not balanced in the region, it implies that part of the flux 
distribution is missing from the limited-size of the magnetogram or that the 
cahbration of the magnetograph is bad. In both cases, any magnetic extrap- 
olation method will only derive an approximate field. Unfortunately the size 
of observed vectormagnetograms is limited and a real magnetogram wiU 
usually not be exactly flux balanced. Field lines passing the boundary of the 
computational box can be either open field lines or close on the photosphere 
outside the observed magnetogram. 
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not necessarily vanish. It is straightforward to extend the 
iteration by 

dB 

^=mG (36) 

on the open boundaries. Equation (36) changes the bound- 
ary values in such way that L decreases. In the following we 
choose the lateral and top boundaries as open (they are ini- 
tiaUzed with a potential field) and apply Eq. (36) here. The 
bottom boundary remains time independent during the iter- 
ation. Let us remark that one might as well use the above 
described photospheric boundary condition for isolated ac- 
tive regions (j„ fixed on the photosphere only in regions of 
positive Bn). To do so one has to apply Eq. (36) also for the 
transversal magnetic field on the photosphere in regions with 
negative Bn- 

If we do not keep the lateral and top boundary conditions 
fixed, but allow them to relax with help of Eq. (36), L drops 
to 36.2 which is only one order of magnitude above the dis- 
cretisation error. The corresponding value of the total force 
decreases sUghtly (|j x B| = 5.3) compared to fixed bound- 
ary conditions. For observed vectormagnetograms it might 
be useful to choose a sufficiently large area around an active 
region, i.e. with the side boundaries relatively far away from 
the non-potential active region, so that the magnetic field can 
be assumed to be approximately potential close to the side 
boundaries. We remark that it might be hard to find cor- 
responding vector magnetograms as these instruments (e.g. 
IVM in Hawaii) do not observe the complete solar disk, but 
only restricted areas. We have also carried out a run on the 
80 X 80 X 40 grid with potential field boundary conditions. 
For this run a stationary state is reached after 10, 000 steps, 
and the forces are less than one order of magnitude higher 
than the discretisation error, while L is nearly three orders of 
magnitude above the discretisation error. Please note that the 
absolute error in the |j x B | has the same order of magnitude 
here as for the corresponding 40 x 40 x 20 run. If we allow 
relaxation on the side and top boundaries, L drops to 96 and 
|j X B| = 8.7 is slightly higher than for fixed side boundary 
conditions. We conclude that the uncertainty of the side and 
top boundary conditions affects the convergence of the opti- 
mization method. Its influence can be reduced significantly if 
we allow for a relaxation of these boundaries and only fix the 
photospheric boundary conditions, but further improvements 
would still be welcome. 

5 Conclusions 

In the present paper we have presented an assessment of 
the properties of an optimization method which has recently 
been proposed for the calculation of nonhnear force-free fields 
from given boundary data (Wheatland et al. , 2000). One part 
of the assessment was a comparison the performance of the 
optimization method with the performance of an MHD relax- 
ation method. For both methods new parallelized codes have 
been developed. Both methods have been applied to finding 



a known semi-analytic solution (Low and Lou , 1990) from a 
given non-equilibrium initial condition. Both methods con- 
verge to the exact solution, but the optimization method has 
higher accuracy. The MHD relaxation method is only applied 
to this case since in its present numerical implementation, the 
boundary conditions could give rise to inconsistencies. 

To simulate a more reahstic situation which is closer to 
working with observational data from vector magnetograms 
we added noise to the boundary conditions. We have only 
applied the optimization method to this type of problem and 
found that already a relatively small noise level of 10 % 
can affect the convergence of the method quite considerably. 
More work is needed to see how this difficulty can be over- 
come. We also intend to generalize the relaxation method in 
such way that is able of coping with more realistic boundary 
data. 

The ultimate aim for the future is to apply these methods to 
data from vector magnetograms. Unfortunately vector mag- 
netograms are less accurate than line-of-sight magnetograms 
(e.g. from MDI on SOHO). Therefore, in the hght of our re- 
sults in the cases where noise was added to the boundary con- 
ditions it will be interesting to see whether and how quickly 
the methods converge. We also plan to use stereoscopic in- 
formation as a further constraint in the reconstruction pro- 
cess. This will become especially important for the analy- 
sis of data from the STEREO mission. A method which is 
based on linear force free fields has recently been proposed 
by Wiegelmann and Neukirch (2002), but hnear force free 
fields seem to be too restrictive to describe coronal phenom- 
ena appropriately. Therefore, the natural next step will be the 
use of nonhnear force free fields in such a code. 
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